clearvars
clc

Results      = NaN(25,9);
Results(:,1) = linspace(2001,2025,25)';

for ical=1:2

    name = strcat('Results',num2str(ical));
    load(name)

    % Construct some required objects
    
    LDCfa                = LHCfa .* LDBas;
    LLBas                = lLBas;
    LLBas(:,2:end,2:end) = repmat(lLBas(:,2:end,1),1,1,T-1).*cumprod(LDBas,3);
    LLCfa                = lLCfa;
    LLCfa(:,2:end,2:end) = repmat(lLCfa(:,2:end,1),1,1,T-1).*cumprod(LDCfa,3);
    bartikjose           = xlsread('Inputs/InputData.xlsx','EXP','B2:B51');
    barr                 = bartikjose*2.63/mean(bartikjose);
    
    % Define baseline quantities
    
    PopuB              = squeeze(sum(lLBas(1:M,:,:),2));
    TotalLaborSupplyB  = squeeze(sum(lLBas(1:M,2:end,:),2));
    TotalEmploymentB   = squeeze(sum(LLBas(1:M,2:end,:),2));
    MoPB               = squeeze(sum(LLBas(1:M,2:13,:),2))./PopuB;
    NoPB               = squeeze(sum(LLBas(1:M,14:15,:),2))./PopuB;
    UoPB               = (TotalLaborSupplyB-TotalEmploymentB)./PopuB;
    loPB               = TotalLaborSupplyB./PopuB;
    
    % Define counterfactual quantities
    
    PopuC              = squeeze(sum(lLCfa(1:M,:,:),2));
    TotalLaborSupplyC  = squeeze(sum(lLCfa(1:M,2:end,:),2));
    TotalEmploymentC   = squeeze(sum(LLCfa(1:M,2:end,:),2));
    MoPC               = squeeze(sum(LLCfa(1:M,2:13,:),2))./PopuC;
    NoPC               = squeeze(sum(LLCfa(1:M,14:15,:),2))./PopuC;
    UoPC               = (TotalLaborSupplyC-TotalEmploymentC)./PopuC;
    loPC               = TotalLaborSupplyC./PopuC;
    
    % Run ADH regressions
    
    betabig            = zeros(4,24);
    
    for tin=2:T-1
        MoPchan            = (mean(MoPC(:,tin-1:tin+1),2)-mean(MoPB(:,tin-1:tin+1),2))*100;
        NoPchan            = (mean(NoPC(:,tin-1:tin+1),2)-mean(NoPB(:,tin-1:tin+1),2))*100;
        UoPchan            = (mean(UoPC(:,tin-1:tin+1),2)-mean(UoPB(:,tin-1:tin+1),2))*100;
        loPchan            = (mean(loPC(:,tin-1:tin+1),2)-mean(loPB(:,tin-1:tin+1),2))*100;
        X                  = [ones(M,1),barr];
        Y                  = [MoPchan,NoPchan,UoPchan,-loPchan]*(10/(tin-1));
        coeff              = (X'*X)^(-1)*(X'*Y);
        betabig(:,tin)     = coeff(2,:)';
    end
    
    if caltype==1
        Results(:,6:9)=betabig(:,2:26)';
    else
        Results(:,2:5)=betabig(:,2:26)';
    end

end

disp(compose('%.7f', Results))

figure(1)
subplot(2,2,1)
axis([2006 2020 -Inf Inf])
hold on
plot(Results(:,1),Results(:,2),'--diamond','Color','green','LineWidth',1.75)
plot(Results(:,1),Results(:,6),'-^','Color','black','LineWidth',1.75)
xline(2007,'r--','LineWidth',1.75)
hold off
title('\fontsize{12} Manuf. employment/Pop')
xlabel('\fontsize{12} Year')
legend('Model Shock 2011','Model Shock 2007','Location','SouthEast')
subplot(2,2,2)
axis([2006 2020 -Inf Inf])
hold on
plot(Results(:,1),Results(:,3),'--diamond','Color','green','LineWidth',1.75)
plot(Results(:,1),Results(:,7),'-^','Color','black','LineWidth',1.75)
xline(2007,'r--','LineWidth',1.75)
hold off
title('\fontsize{12} Non-manuf. employment/Pop')
xlabel('\fontsize{12} Year')
legend('Model Shock 2011','Model Shock 2007','Location','SouthEast')
subplot(2,2,3)
axis([2006 2020 -Inf Inf])
hold on
plot(Results(:,1),Results(:,5),'--diamond','Color','green','LineWidth',1.75)
plot(Results(:,1),Results(:,9),'-^','Color','black','LineWidth',1.75)
xline(2007,'r--','LineWidth',1.75)
hold off
title('\fontsize{12} NILF/Pop')
xlabel('\fontsize{12} Year')
legend('Model Shock 2011','Model Shock 2007','Location','NorthEast')
subplot(2,2,4)
axis([2006 2020 -Inf Inf])
hold on
plot(Results(:,1),Results(:,4),'--diamond','Color','green','LineWidth',1.75)
plot(Results(:,1),Results(:,8),'-^','Color','black','LineWidth',1.75)
xline(2007,'r--','LineWidth',1.75)
hold off
title('\fontsize{12} Unemployment/Pop')
xlabel('\fontsize{12} Year')
legend('Model Shock 2011','Model Shock 2007','Location','NorthEast')
